## RSID: The interpolating decoder for RS codes

from PyM import *

def primitive_root(K):
    q = cardinal(K); n = q-1
    for j in range(2,q):
        if order(element(j,K))==n:
            return element(j,K)
    return "primitive_root error: root not found"
    

def RSID(a,k,y):
    n = len(a); t = (n-k)//2
    P = vandermonde(a,n-t); Q = vandermonde(a,n-k-t+1)
    for j in range(n):
        Q[:,j] = y[j] * Q[:,j] 
    PQ = left_kernel(stack(P,Q))[0]
    P = PQ[:(n-t)]; Q = PQ[(n-t):]
    if is_zero(Q): return "RSID Error: Q = 0"
    [_,T] = polynomial_ring(K_(a),'T')
    P = polynomial(P,T); Q = polynomial(Q,T)
    if not P % Q: return "RSID Error: Q does not divide P"
    f = -P//Q
    return vector([evaluate(f,t) for t in a])

q=23   # Under 1 s
#q=101 # Over 20 s
K = Zn(q)
w = primitive_root(K)
n=q-1
a = geometric_series(w,n)

k = (3*n)//4
if k%2: k=k+1

t = (n-k)//2
e = rd_error_vector(n,t,K)

show("error vector =",e)


x = RSID(a,k,e)
show(x)